home *** CD-ROM | disk | FTP | other *** search
/ Amiga Packmags / Source, The - Issue 5 (1993)(Epsilon)[WB].zip / Source, The - Issue 5 (1993)(Epsilon)[WB].adf / Source / Vectors / VectorSpace / PointInPoly.lha / pointnpolygon.txt3 < prev   
Internet Message Format  |  1989-10-24  |  6KB
  1. From albanycs!leah:rsb584 Tue Mar  1 15:54:24 1988
  2. Received: by albanycs.albany.edu (5.54/4.8)
  3.     id AA11798; Tue, 1 Mar 88 15:48:08 EST
  4. Date: Tue, 1 Mar 88 15:48:05 EST
  5. From: albanycs!leah:rsb584 (Raymond S Brand)
  6. Received: by leah.Albany.EDU (5.58/1.1)
  7.     id AA15407; Tue, 1 Mar 88 15:48:05 EST
  8. Message-Id: <8803012048.AA15407@leah.Albany.EDU>
  9. To: albanycs:beowulf!rsbx
  10. Subject: winding.num.txt
  11.  
  12. >From kenm@sci.UUCP Sat Feb 27 19:21:05 1988
  13. Path: leah!uwmcsd1!tut.cis.ohio-state.edu!mailrus!nrl-cmf!ames!oliveb!pyramid!voder!tolerant!sci!kenm
  14. From: kenm@sci.UUCP (Ken McElvain)
  15. Newsgroups: comp.graphics,comp.sys.ibm.pc
  16. Subject: Re: point-in-polygon again
  17. Summary: point in poly by winding number
  18. Keywords: READ THIS ONLY IF YOUR IQ <<125
  19. Message-ID: <16460@sci.UUCP>
  20. Date: 28 Feb 88 00:21:05 GMT
  21. References: <971@ut-emx.UUCP> <20533@amdcad.AMD.COM> <7626@pur-ee.UUCP> <5305@spool.cs.wisc.edu>
  22. Distribution: na
  23. Organization: Silicon Compilers Systems Corp. San Jose, Ca
  24. Lines: 132
  25.  
  26. In article <5305@spool.cs.wisc.edu>, planting@speedy.cs.wisc.edu (W. Harry Plantinga) writes:
  27. > In article <7626@pur-ee.UUCP> beatyr@pur-ee.UUCP (Robert Beaty) writes:
  28. > >In article <20533@amdcad.AMD.COM> msprick@amdcad.UUCP (Rick Dorris) writes:
  29. > >>In article <971@ut-emx.UUCP> jezebel@ut-emx.UUCP (Jim @ MAC) writes:
  30. > >>>Have a nice one here:
  31. > >>>Have a boundary defined on the screen. Boundary is composed of points
  32. > >>>joined by lines... Now some random points are generated and I want to check
  33. > >>>if a given point is within or outside the existing boundary.. Any algorithm for
  34. > >> <rick suggests the infinite ray/intersection counting algorithm>
  35. > >
  36. > >I have seen this type of algorithm before and the one I thought up
  37. > >in an afternoon is FAR surperior and vastly simpler to code up.
  38. > >
  39. > ><robert beaty suggests an algorithm that measures sums the angles of
  40. > >the lines connecting the poitns and the test point>
  41. > 
  42. > Haven't I seen this discussion of point-in-polygon algorithms about 15
  43. > times before? :-)
  44. > 
  45. > Robert, your algorithm is simpler only in the kind of problems it
  46. > handles.  It will only work for convex polygons, whereas the other
  47. > works for any polygons.  It's not even much simpler; measuring angles
  48. > and determining whether line segments intersect are of comparable
  49. > difficulty.
  50. > 
  51. > Maybe you should have said that althogh your algorithm is far
  52. > inferior, it's not any easier to code?
  53.  
  54.  
  55. Robert's suggestion is a good one.  The sum of the angles about the
  56. test point is known as the winding number.  For non intersecting
  57. polygons if the winding number is non-zero then the test point is
  58. inside the polygon.  It works just fine for convex and concave
  59. polygon's.  Intersecting polygon's give reasonable results too.
  60. A figure 8  will give a negitive winding number for a test point
  61. in one of the loops and a positive winding number for the other loop,
  62. with all points outside having a winding number of 0.  Other advantages
  63. of the winding number calculation are that it does not suffer from
  64. the confusion of the infinite ray algorithm when the ray crosses a vertex
  65. or is colinear with an edge.
  66.  
  67. Here is my version of a point in poly routine using a quadrant
  68. granularity winding number.  The basic idea is to total the angle
  69. changes for a wiper arm with its origin at the point and whos end
  70. follows the polygon points.  If the angle change is 0 then you are
  71. outside, otherwise you are in some sense inside.  It is not necessary to
  72. compute exact angles, resolution to a quadrant is sufficient, greatly
  73. simplifying the calculations. 
  74.  
  75. I pulled this out of some other code and hopefully didn't break it in
  76. doing so.  There is some ambiguity in this version as to whether a point
  77. lying on the polygon is inside or out.  This can be fairly easily
  78. detected though, so you can do what you want in that case. 
  79.  
  80.  
  81. -----------------------------------------------------------------
  82. /*
  83.  * Quadrants:
  84.  *    1 | 0
  85.  *    -----
  86.  *    2 | 3
  87.  */
  88.  
  89. typedef struct  {
  90.         int x,y;
  91. } point;
  92.  
  93. pointinpoly(pt,cnt,polypts)
  94. point pt;       /* point to check */
  95. int cnt;        /* number of points in poly */
  96. point *polypts; /* array of points, */
  97.                 /* last edge from polypts[cnt-1] to polypts[0] */
  98. {
  99.         int oldquad,newquad;
  100.         point thispt,lastpt;
  101.         int a,b;
  102.         int wind;       /* current winding number */
  103.  
  104.         wind = 0;
  105.         lastpt = polypts[cnt-1];
  106.         oldquad = whichquad(lastpt,pt); /* get starting angle */
  107.         for(i=0;i<cnt;i++) { /* for each point in the polygon */
  108.                 thispt = polypts[i];
  109.                 newquad = whichquad(thispt,pt);
  110.                 if(oldquad!=newquad) { /* adjust wind */
  111.                         /*
  112.                          * use mod 4 comparsions to see if we have
  113.                          * advanced or backed up one quadrant 
  114.                          */
  115.                         if(((oldquad+1)&3)==newquad) wind++;
  116.                         else if((((newquad+1)&3)==oldquad) wind--;
  117.                         else {
  118.                                 /*
  119.                                  * upper left to lower right, or
  120.                                  * upper right to lower left. Determine
  121.                                  * direction of winding  by intersection
  122.                                  *  with x==0.
  123.                                  */
  124.                                 a = lastpt.y - thispt.y;
  125.                                 a *= (pt.x - lastpt.x);
  126.                                 b = lastpt.x - thispt.x;
  127.                                 a += lastpt.y * b;
  128.                                 b *= pt.y;
  129.  
  130.                                 if(a > b) wind += 2;
  131.                                 else wind -= 2;
  132.                         }
  133.                 }
  134.                 lastpt = thispt;
  135.         }
  136.         return(wind); /* non zero means point in poly */
  137. }
  138.  
  139. /*
  140.  * Figure out which quadrent pt is in with respect to orig
  141.  */
  142. int whichquad(pt,orig)
  143. point pt;
  144. point orig;
  145. {
  146.         if(pt.x < orig.x) {
  147.                 if(pt.y < orig.y) quad = 2;
  148.                 else quad = 1;
  149.         } else {
  150.                 if(pt.y < orig.y) quad = 3;
  151.                 else quad = 0;
  152.         }
  153. }
  154.  
  155.  
  156. Ken McElvain
  157. decwrl!sci!kenm
  158.  
  159.  
  160.